
phi_id0  = para(1);
phi_id1  = para(2);
phi_sig0 = para(3);
phi_sig1 = para(4);
phi_mu   = para(5);
a        = para(6);
b        = para(7);

gamma   = para(8);
psi     = 1.5;

rho_CE  = 0.2; 

p       = 0.878776;

mu1     = 0.003060;
mu2     = 0.000453;

sig1    = 0.001500;
sig2    = 0.002768;
sig3    = 0.007034;

DD           = 44; % mean disaster duration in month
i            = 1:100000;
duration_fun = @(p)sum((1-p)*p.^(i-1).*i);
p3          = fzero(@(p)duration_fun(p)-DD,0.99);

ED          = 0.834185141; % mean drop in consumption (simple growth rate)
i           = 1:100000;
mean_fun    = @(m3)sum((1-p3)*p3.^(i-1).*exp(i*m3+i*sig3^2/2));
mu3         = fzero(@(m)mean_fun(m)-ED,-0.0038);

pd          = fzero(@(x)solve_P(x,x,p,p,p3)-(44/12)/(2019-1929),1e-3); % This is pi in the paper.

mu_C        = [mu1;mu2;mu3];
sig_C       = [sig1;sig2;sig3];

P = [ p*(1-pd)      (1-p)*(1-pd)  pd;
      (1-p)*(1-pd)  p*(1-pd)      pd;
      (1-p3)/2      (1-p3)/2      p3];

P_nodis = [p,1-p;1-p,p];

nS      = 3;

% ergodic distribution
stat        = P^1e6;
stat        = stat(1,:)';


% Computation
Tmax        = 12;   % maximum CDS maturity (in months)
Tmax_CDS    = 60;   % maximum CDS maturity (in months)
Tmax_IV     = 12;
k_min       = 0.1;
k_max       = 200;
nK          = 8000;
nK_fine     = 1*nK;
nE          = 51;

meth        = 'linear';
ex_meth     = 'linear';

% Convergence for VFI
prec        = 1e-6;


% Frictions
tau_i       = 0.296; % tax rate on interest income
tau_d       = 0.112; % tax rate on dividend income
tau_c       = 0.329; % tax rate on corporate earnings

xi_d        = 0.01; % debt issuance fee
xi_e        = 0.04; % equity issuance costs

sig_Xid     = phi_id0  + phi_id1  * (sig_C-stat'*sig_C);
sig_XS      = phi_sig0 + phi_sig1 * (sig_C-stat'*sig_C);

sig_X       = sqrt( sig_XS.^2 + sig_Xid.^2 );

mu_X        = (mu_C-stat'*mu_C)*phi_mu + stat'*mu_C;
mu_XQ       = mu_X - gamma*rho_CE*sig_XS.*sig_C;

SR          = mu_C./sig_C;
E_SR        = SR'*stat;
SR          = SR - E_SR;

omega       = a - b*SR;
omega       = min(1,max(omega,0));

